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Abstract 

This  document  is  the  final  report  of  the  research  sponsored  by  AFOSR  under  contract  No. 
F49620-91-C-0048-P00002  with  Stanford  University.  This  research  effort  focused  on  the  math¬ 
ematical  modeling  and  numerical  simulation  of  the  dynamics  of  flexible  structures  undergoing 
large  overall  motions.  The  results  reported  below  fall  within  the  following  areas:  (i)  Enhanced 
numerical  models  for  shell  intersections,  (ii)  Asymptotic  methods  for  nonlinear  shells,  and  (iii) 
Conserving  time-stepping  algorithms  for  nonlinear  dynamics. 

1.  Introduction 

Realistic  models  of  spacecrafts,  space  platforms  and  related  devices  invariably  involve 
large  and  very  flexible  structures  undergoing  large  overall  motions,  thus  rendering  inappro¬ 
priate  analyses  based  on  a  linearized  approach.  The  ultimate  goal  of  this  investigation  is 
the  development  of  analytical  and  computational  techniques  for  the  mathematical  model¬ 
ing  and  numerical  simulation  of  the  nonlinear  response  of  this  type  of  structure  undergoing 
such  large  overall  motions.  Central  to  the  proposed  approach  is  the  use  of  nonlinear  me¬ 
chanical  models  which  place  no  restriction  on  the  allowable  flexibility  of  the  structure, 
and  capture  dynamic  coupling  effects  induced  by  the  deformation  of  the  structure  without 
resorting  to  ad-hoc  simplifications.  In  this  way,  a  crucial  feature  of  the  mathematical 
models  under  consideration  is  the  rich  geometric  structure  associated  to  them.  The  results 
obtained  in  this  research  project  show  that  the  exploitation  of  this  geometric  structure 
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leads  to  very  efficient  numerical  techniques  which  exhibit  an  excellent  computational  per¬ 
formance. 

The  research  developed  under  this  project  has  been  the  continuation  of  the  research 
sponsored  by  AFOSR  under  grants  No.  DJA/AFOSR  86-0292  and  No.  2DJA/AFOSR  89- 
0294  with  Stanford  University  during  the  period  1986-1990.  Accomplished  results  of  these 
previous  projects  were:  (i)  The  formulation  and  analysis  of  geometrically  exact  models, 
(ii)  The  numerical  analysis  of  the  dynamics  of  shell  models,  and  (iii)  The  formulation 
of  the  Energy-Momentum  method  for  the  stability  analysis  of  Hamiltonian  systems.  As 
indicated  in  the  original  proposal,  the  main  goals  of  the  present  project  are  the  extension 
and  generalization  of  these  results  according  to  the  following  lines  of  research: 

i.  Formulation  and  numerical  analysis  of  enhanced  models  for  junctions  and  shell  in¬ 
tersections  in  the  framework  of  the  previously  developed  geometrically  exact  models. 

ii.  Asymptotic  methods  for  nonlinear  shells. 

iii.  Formulation  of  conserving  time-stepping  algorithms  for  nonlinear  dynamics  in  the 
context  of  nonlinear  finite  element  methods. 

The  results  obtained  along  these  three  lines  of  research  are  summarized  in  the  following 
section.  See  Simo  [1992],  Fox,  Raoult  &  Simo  [1993],  and  Simo,  Rifai  &  Fox  [1992] 
for  a  complete  description  of  these  results  in  the  above  three  respective  areas. 


2.  Summary  of  the  Main  Accomplishments 

This  section  presents  the  main  results  obtained  in  this  research  project.  These  ac¬ 
complishments  are  grouped  in  the  three  main  lines  of  research  indicated  above. 

2.1.  Enhanced  models  for  junction  and  shell  intersections 

As  noted  in  the  introduction,  previous  research  efforts  under  the  support  of  the 
AFOSR  concentrated  in  the  formulation  and  numerical  analysis  of  geometrically  nonlinear 
shell  models.  See  SlMO  &  Fox  [1989],  SlMO,  Fox  &  Rifai  [1989,1990]  and  SlMO,  Fox 
&  RlFAI  [1990]  for  a  complete  description.  One  of  the  original  goals  accomplished  in 
the  present  research  effort  has  been  the  extension  and  generalization  of  these  models  to 
accommodate  situations  arising  in  the  actual  design  of  space  structures.  In  particular, 
non-smooth  intersections  of  shells  appear  almost  invariably  in  such  designs.  The  complete 
mathematical  characterization  and  numerical  simulation  of  their  mechanical  response  is 
then  of  major  practical  interest. 

In  recent  years,  there  has  been  a  significant  activity  in  computational  treatments  of 
nonlinear  shell  theory.  Representative  examples  are  the  recent  works  of  PARISH  [1991], 
where  the  classical  degenerated  solid  concept  is  revisited  in  a  nonlinear  setting  that  in¬ 
corporates  geometrically  exact  rotational  updates,  BUCHTER  &  Ramm  [1992]  where  the 
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FIGURE  2.1.1  Cylindrical  bending  of  a  right  angle  shell,  with  thickness  0.02  units. 
Shown  in  the  figure  are  different  deformed  shapes.  Perfect  agreement  is  found  with 
the  elementary  beam  solution. 


relations  between  classical  shell  theories  (of  the  type  considered  in  this  work)  and  the  de¬ 
generated  solid  approach  are  re-examined  and  further  explored,  and  SANSOUR  &  BUFLER 
[1992]  where  nonlinear  shell  theories  based  on  the  Biot  formulation  introduced  in  ATLURI 
[1984]  are  explored  in  a  computational  context.  An  overview  of  representative  approaches 
in  computational  shell  theory  up  to  the  end  of  the  80’s,  is  contained  in  Noor,  BELYTSCHKO 
&  SlMO  [1990].  Inspection  of  the  large  body  of  literature  concerned  with  computational 
shell  analysis  reveals  that  most  of  the  approaches  are  restricted  by  the  implicit  assumption 
that  the  shell  mid-surface  is  to  remain  smooth.  The  same  limitation  typically  applies  to 
standard  expositions  of  nonlinear  shell  theory;  see  e.g.,  NAGHDI  [1972].  As  a  result  of  this 
restriction,  a  unique  director  field  is  attached  to  each  point  of  the  mid-surface,  which  is 
parametrized  by  two  independent  (rotational)  coordinates. 

It  has  been  widely  accepted  for  a  number  of  years  now  that  the  treatment  of  shell 
intersections  using  finite  elements  requires  formulations  possessing  6  DOF/Node.  Current 
finite  element  methods  of  this  type  incorporate  the  so-called  drill-rotation  about  the  direc¬ 
tor  field  and  often  require  stabilization  via  the  addition  of  an  artificial  drill-stiffness.  The 
recent  work  of  PARISH  [1991]  is  representative  of  this  approach.  Finite  element  methods 
designed  to  incorporate  this  additional  rotation  in  the  membrane  field  have  been  discussed 
in  Allman  [1984],  Taylor  &  Simo  [1985],  Hughes  &  Brezzi  [1989],  Ibrahimbegovic, 
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Taylor  &  Wilson  [1990]  and  Simo,  Fox  &  Hughes  [1992]  among  others.  It  appears 
that  non-classical  theories  which  incorporate  an  independent  rotation  field,  of  the  type 
considered  by  TOUPIN  [1964],  REISSNER  [1965]  and  others,  provide  the  proper  continuum 
setting  for  the  incorporation  of  drill-rotations;  see  e.g.,  Fox  &  SlMO  [1992]. 

In  contrast  with  the  aforementioned  works,  a  direct  approach  based  on  classical  shell 
theory  has  been  proposed  in  this  work  that  circumvents  the  need  for  drill  rotations.  The 
key  observation  exploited  in  the  continuum  formulation  of  the  shell  intersection  problem 
is  the  view  of  a  shell  intersection  as  a  smooth  curve  in  the  reference  shell  mid- surface 
across  which  the  director  field  experiences  a  jump  discontinuity.  Accordingly,  two  (or 
more)  directors  are  assigned  to  each  point  of  an  intersection  in  the  reference  configuration. 
The  additional  assumption  of  a  rigid  shell  intersection  translates  into  the  constraint  which 
requires  that  the  angle  between  the  directors  at  each  point  of  the  intersection  remain 
constant.  Satisfaction  of  this  constraint  can  be  achieved  only  by  demanding  continuity  of 
the  angular  velocity  of  the  director  field  across  the  shell  intersection  (but  not  continuity 
of  the  director  velocity  itself).  Since  the  director  angular  velocity  at  a  point  in  a  shell 
intersection  generally  possesses  three  independent  components,  as  opposed  to  only  two 
independent  components  in  the  case  of  a  smooth  shell  mid-surface,  one  necessarily  obtains 
finite  element  discretizations  possessing  6  global  DOF  at  points  in  a  shell  interaction.  This 
added  degree  of  freedom,  however,  is  totally  unrelated  to  the  drill-rotation  widely  used  to 
tackle  the  problem,  as  indicated  above. 

In  this  context,  a  convenient  parametrization  is  defined  in  terms  of  a  rotation  matrix 
assigned  to  each  node  and  updated  iteratively.  Exact  nonlinear  geometric  updates  are 
employed  for  the  rotational  degrees  of  freedom.  This  procedure  leads  to  a  simple  modifi¬ 
cation  of  previously  developed  finite  element  methods.  In  particular,  at  the  element  level 
the  finite  element  method  remains  unchanged  from  the  standard  formulations  employing 
5  DOF /Node.  The  element  contributions  for  the  residual  vector  and  tangent  matrix  are 
then  assembled  in  the  corresponding  global  arrays,  where  the  different  number  of  degrees 
of  freedom  for  the  shell  intersections  are  taken  in  account  by  the  rotation  matrix  indicated 
above.  This  leads  to  a  simple  and  effective  computational  procedure  that  has  exhibited  an 
excellent  performance  in  a  number  of  numerical  simulations.  Figure  2.2.1  depicts  a  rep¬ 
resentative  numerical  simulation  obtained  with  the  proposed  approach  consisting  of  the 
(nonlinear)  rolling  of  a  L-shape  shell;  see  SlMO  [1992]  for  further  details. 


2.2.  Asymptotic  methods  for  nonlinear  shells 

Plates  and  shells  are  three-dimensional  continuum  bodies  characterized  by  a  reference 
placement  possessing  one  dimension,  the  thickness,  which  is  small  in  comparison  with  the 
other  two  dimensions.  This  distinctive  feature  suggests  the  use  of  an  asymptotic  analysis, 
with  the  thickness  regarded  as  a  small  parameter,  as  a  means  of  justifying  two-dimensional 
mathematical  models  of  plates  and  shells  from  the  full  three-dimensional  theory.  A  large 
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body  of  literature  concerned  with  the  derivation  of  classical  plate  and  shell  equations  does 
in  fact  adopt  asymptotic  expansions' as  a  central  tool  in  the  analysis.  Most  of  the  early 
derivations  or  justifications  of  plate  and  shell  models  via  asymptotic  analysis  are  typically 
restricted  to  the  linear  theory. 

A  reformulation  and  precise  formalization  of  the  asymptotic  approach  within  the 
general  framework  of  three-dimensional  nonlinear  elasticity  is  described  in  ClARLET  & 
DESTUYNDER  [1979]  and  ClARLET  [1980].  There,  it  is  shown  that  the  von  Karman  model 
of  plates  arises  as  the  leading  term  in  a  formal  asymptotic  expansion  of  the  nonlinear 
three-dimensional  theory  cast  as  a  mixed  variational  problem  in  stresses  and  displace¬ 
ments.  An  alternative  formulation  of  this  approach,  which  uses  a  variational  formulation 
solely  in  terms  of  displacements,  was  given  by  RAOULT  [1988].  The  asymptotic  approach 
has  been  applied  by  these  authors  and  their  co-workers  to  a  variety  of  models,  including 
shallow  shells  and  junctions.  The  lack  of  invariance  under  superposed  rigid  body  motions 
(Euclidean  isometries)  exhibited  by  the  (second  order)  models  arising  in  this  asymptotic 
approach,  the  von  Karman  model  in  particular,  is  a  noteworthy  feature  also  shared  by  the 
linear  theory  of  elasticity.  This  striking  feature  is  to  be  contrasted  with  the  full  invariance 
properties  under  Euclidean  isometries  of  the  nonlinear  three-dimensional  theory  which  is 
taken  as  the  point  of  departure  in  these  asymptotic  analyses. 

An  entirely  different  approach  to  the  formulation  of  one  and  two-dimensional  math¬ 
ematical  models  for  rods,  plates  and  shells  employs  directly  the  notion  of  directed  media 
introduced  in  COSSERAT  [1909].  This  approach  was  revitalized  in  the  pioneering  paper 
of  ERICKSEN  &  TRUESDELL  [1958]  and  encompasses  classical  rod  models  going  back  to 
EULER,  Clebch  &  KlRCHHOFF  (see  e.g.,  Love  [1932]).  Comprehensive  reviews  of  this 
methodology,  often  referred  to  as  the  direct  approach ,  can  be  found  in  Antman  [1972] 
and  NAGHDI  [1972]. .  Two  representative  examples  are  the  constrained  two-director  rod 
model,  which  incorporates  finite  extension,  shearing,  bending  and  twist,  and  the  inexten- 
sible  one-director  shell  model,  which  incorporates  finite  membrane  and  bending  as  well 
as  transverse  shear  deformation.  As  stressed  in  ANTMAN  [1976],  a  significant  and  often 
overlooked  property  of  models  arising  in  projection-constraint  methods  (intimately  related 
to  the  direct  approach)  is  the  exact  character  of  the  balance  laws  of  momentum,  which  can 
be  derived  directly  from  the  three-dimensional  theory.  Furthermore,  these  models  enjoy 
the  same  invariance  properties  under  superposed  Euclidean  isometries  as  the  full  nonlinear 
three-dimensional  theory. 

From  the  above  remarks,  one  would  be  tempted  to  conclude  that  the  von  Karman 
model  of  plates,  and  not  any  of  the  full  nonlinear  models  obtained  by  a  direct  approach,  is 
the  canonical  nonlinear  model  that  arises  as  the  leading  term  in  an  asymptotic  expansion 
of  the  three-dimensional  theory.  It  was  shown  in  this  work  that  full  nonlinear  and  properly 
invariant  models  of  plates,  referred  to  as  the  membrane  theory  and  the  inextensional  theory, 
do  in  fact  arise  as  leading  terms  in  an  asymptotic  expansion  of  the  three-dimensional 
nonlinear  theory.  Our  approach  relied  on  a  reformulation  of  the  asymptotic  procedure 
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introduced  in  ClARLET  [1980]  and  further  developed  in  RAOULT  [1988],  cast  in  a  form 
designed  to  preserve  the  full  invariance  properties  of  the  three-dimensional  theory.  Only 
two  assumptions  were  made  in  our  formal  asymptotic  analysis: 

i.  The  applied  loads  admit  an  expansion  in  power  series  of  the  thickness  £,  after  a 
reparametrization  to  an  £-independent  domain  is  introduced. 

ii.  The  ansatz  is  made  that  the  deformations  of  the  plate,  reparametrized  to  the 
^-independent  domain,  can  be  formally  expanded  in  powers  of  e. 

We  showed  that,  within  the  hierarchy  of  nonlinear  plate  models  arising  in  the  asymp¬ 
totic  procedure,  the  appearance  of  a  specific  model  as  the  leading  term  is  dictated  only 
by  the  order  of  the  leading  non-zero  term  in  the  expansion  of  the  loads.  In  particular, 
we  showed  that  the  nonlinear  membrane  theory  arises  if  the  body  forces  and  edge  loading 
are  order  1  and  the  surface  loading  is  order  e.  In  this  classical  model,  the  stored  energy 
function  depends  only  on  the  first  fundamental  form  of  the  deformed  plate  mid- surface 
and  the  solutions  are  then  characterized  as  minimizers  of  the  total  potential  energy.  No 
restrictions  are  placed  by  the  model  on  the  magnitude  of  the  membrane  strains. 

Further,  we  showed  that  the  nonlinear  inextensional  theory  arises  as  the  leading  term 
in  the  formal  asymptotic  expansion  for  body  forces  and  edge  loading  of  order  £,  surface 
loading  of  order  e 2  and  vanishing  transverse  normal  force  (i.e.,  pure  bending  loads).  In  this 
model,  the  stored  energy  function  depends  only  on  the  second  fundamental  form  of  the 
deformed  plate  mid-surface,  which  must  be  isometric  to  a  plane  since  the  first  fundamental 
form  is  the  identity.  The  solutions  of  the  problem  are  characterized  via  the  asymptotic 
procedure  as  minimizers  of  the  total  potential  energy  over  the  manifold  of  surfaces  isometric 
to  the  plane.  No  restrictions  are  placed  on  the  magnitude  of  the  bending  strains.  Again 
it  is  shown  that  precise  restrictions  on  the  order  of  magnitude  of  the  stress  field  are  not 
postulated  a  priori  (as  in  previous  works)  but  are  obtained  as  a  result  in  the  course  of  the 
asymptotic  derivation. 

Finally,  we  showed  that  the  von  Karman  model  of  plates  arises  as  the  leading  term 
in  the  asymptotic  expansion  for  order  e2  in-plane  and  e3  vertical  body  force  and  edge 
loading,  along  with  order  £3  in-plane  and  £4  vertical  surface  loading.  The  preceding  results, 
summarized  in  the  table,  provide  a  precise  description  of  the  position  of  these  theories 
relative  to  the  magnitude  of  the  applied  loads.  In  the  development  of  these  results  the 
simplifying  constitutive  assumption  is  made  of  a  Saint  Venant-KirchhofF  material.  The 
limitations  of  this  model  are  well-known.  We  investigated  the  extension  of  our  results 
to  general  hyperelastic  constitutive  models  and  made  a  preliminary  outline  of  the  steps 
required  to  perform  such  an  analysis. 

The  models  described  above  can  be  derived  from  the  one-director  shell  model  (related 
to  the  direct  approach,  see  SlMO  &  Fox  [1989])  by  placing  suitable  restrictions  on  the 
director  field  and  the  shell  mid-surface.  In  particular,  the  membrane  theory  is  obtained 
in  the  absence  of  a  director  field,  while  the  inextensional  theory  arises  when  the  director 
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TABLE  2.2.1  Scaling  of  the  external  loading  leading  to  different  plate  models. 


Body  and  Edge  Load 

Surface  Load 

Plate  Model 

in-plane 

vertical 

in-plane 

vertical 

Membrane 

1 

1 

£ 

£ 

Inextensional* 

£ 

£ 

£ 2 

£2 

von  Karman 

£3 

e3 

£4 

*  For  the  inextensional  theory,  the  further  restriction  must  be  made  that  no  order  e  resultant  force  exists 
(only  order  e  resultant  moment  loading). 


field  is  identified  with  the  normal  to  the  deformed  mid-surface  which  is  constrained  to 
remain  isometric  to  the  reference  mid-surface  (here,  a  plane).  These  results,  therefore, 
lend  considerable  support  to  the  one-director  shell  model  as  a  general  theory  applicable 
to  nonlinear  plates /shells.  Moreover,  they  provide  a  justification  for  the  success  obtained 
with  this  model  in  direct,  large  scale,  numerical  simulations  of  complex  problems  subjected 
to  general  loading.  For  a  comprehensive  account  of  numerical  analysis  aspects  related  to 
the  implementation  of  the  one-director  nonlinear  shell  model  we  refer  to  SlMO,  RlFAI  & 
Fox  [1990]  and  references  therein. 


2.3.  Conserving  time-stepping  algorithms  for  nonlinear  dynamics 

Of  special  interest  is  the  dynamics  of  very  flexible  shell  structures  undergoing  finite 
deformations  and  large  overall  motions.  In  this  situation  there  are  two  fundamental  con¬ 
stants  of  the  motion,  the  total  linear  and  angular  momentum  maps,  which  are  conserved 
by  the  dynamics  regardless  of  the  specific  functional  form  assumed  for  the  constitutive 
response.  If  one  further  specializes  the  shell  model  by  assuming  hyperelastic  response  and 
conservative  loading,  the  mechanical  system  becomes  Hamiltonian  (see  SlMO,  MaRSDEN 
&  KRISHNAPRASSAD  [1988]),  and  the  total  energy  (i.e.,  the  Hamiltonian)  is  also  conserved 
by  the  dynamics.  A  main  goal  of  this  project  is  the  design  of  time-stepping  algorithms 
which  will  automatically  inherit  these  fundamental  conservation  properties  present  in  the 
continuum  model.  We  showed  in  this  work  that  these  conservation  laws  are  preserved  in 
the  time  discretization  only  by  a  very  specific  class  of  time-stepping  algorithms.  Even 
if  this  class  of  algorithms  is  adopted,  the  conservation  properties  will  not,  in  general,  be 
preserved  by  the  subsequent  spatial  discretization.  A  further  goal  of  this  work  is  the  for¬ 
mulation  of  finite  element  interpolations  which  preserve  the  conservation  laws  inherited  by 
the  time-stepping  algorithm.  Our  main  results  can  be  summarized  as  follows. 

i.  Momentum  conserving  time-stepping  algorithms.  A  crucial  aspect  in  the  numerical 
simulation  of  the  dynamics  of  nonlinear  structural  models  concerns  the  time-intergration 
of  the  rotational  dynamics  which  evolves  in  the  rotation  group.  The  difficulty  is  asso¬ 
ciated  with  the  geometry  of  the  rotation  group  which  is  not  a  (flat)  linear  space  but  a 
differentiable  manifold.  Recently,  SlMO  &  WONG  [1990]  showed  that  the  extension  of  the 
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Newmark  formulae  to  the  rotation  group  written  in  body  coordinates,  in  the  form  sug¬ 
gested  in  SlMO  &  Vu-Quoc  [1988],  along  with  a  generalized  mid-point  rule  approximation 
to  Euler’s  angular  momentum  equation  written  in  conservation  form  (or  impulse  form  in 
the  terminology  of  ZlENKIEWICZ,  WOOD  &  TAYLOR  [1980]),  yields  an  algorithm  which 
exactly  preserves  total  angular  momentum.  In  addition,  the  algorithm  is  frame  invariant 
and  conserves  energy  in  the  absence  of  external  loading  only  for  a  choice  of  parameters 
leading  to  the  mid-point  rule. 

Two  nontrivial  difficulties  arise  in  the  extension  of  the  preceding  class  of  algorithms 
to  the  dynamics  of  nonlinear  shells.  First,  the  coupling  between  the  rotational  and  trans¬ 
lational  parts,  which  manifests  itself  in  the  expression  for  the  total  angular  momentum 
map.  Second,  the  additional  S1  symmetry  present  in  classical  shell  models,  which  is  the 
result  of  the  invariance  of  the  shell  equations  with  respect  to  a  one-parameter  group  of 
(drill)  rotations  about  the  director  field.  As  discussed  in  SlMO  &  Fox  [1989]  and  SlMO, 
Fox  &  RlFAI  [1990],  the  enforcement  of  this  invariance  condition  leads  to  a  replacement  of 
the  rotation  group  by  the  unit  sphere,  since  only  two  independent  parameters  are  required 
to  describe  the  orientation  of  the  director  field.  The  class  of  algorithms  proposed  in  this 
research  is  then  constructed  as  follows: 

1.1.  Both  the  strong  and  weak  forms  of  the  Euler  laws  of  motion  for  nonlinear  shells 
are  written  in  conservation  form  and  approximated  via  a  generalized  mid-point 
rule. 

1.2.  The  Newmark  formulae  are  extended  to  the  unit  sphere  leading  to  a  time¬ 
stepping  procedure  which  exactly  preserves  frame  invariance  and  is  consistent 
with  the  classical  rule  of  composition  of  angular  velocities. 

1.3.  The  configurations  are  updated  using  the  exponential  mapping  in  the  unit 
sphere,  thus  ensuring  exact  satisfaction  of  the  unit  director  constraint  condi¬ 
tion. 

We  showed  that  (for  Neumann  boundary  conditions  and  equilibrated  loading)  total  angular 
momentum  is  exactly  conserved  only  for  one  member  of  the  preceding  family  of  algorithms: 
the  conservation  form  of  the  mid-point  rule.  The  proof  shows  that  this  conservation  prop¬ 
erty  depends  critically  on  adopting  the  steps  i.l  through  i.3  above.  An  identical  result  also 
holds  for  three-dimensional  continuum  mechanics.  In  fact,  for  Newmark-like  algorithms, 
exact  preservation  of  the  conservation  law  of  total  angular  momentum  holds  only  for  a 
mid-point  rule  approximation  to  the  conservation  form  of  the  momentum  equations.  We 
remark  that  the  accelerations  play  no  role  in  such  a  formulation,  and  are  obtained  merely 
via  postprocessing. 

ii.  Energy  conserving  algorithms.  Exact  energy  conserving  algorithms  for  Hamilto¬ 
nian  systems  have  received  considerable  attention  in  the  literature;  see  e.g.  GREENSPAN 
[1974],  Bayliss  &  Isaacson  [1975],  LaBudde  &  Greenspan  [1976a, b]  and  Hughes, 
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LlU  &  CAUGHY  [1978]  among  others.  For  the  linear  problem,  it  is  well-known  that  the 
mid-point  rule  is  the  only  one-step  ‘  method  possessing  this  property;  see  BELYTSCHKO 
&  SCHOEBERLE  [1975].  We  have  proved  a  somewhat  surprising  generalization  of  this  re¬ 
sult:  Given  any  stored-energy  function  and  given  an  arbitrary  potential  for  the  external 
loading,  there  exists  a  generalized  mid-point  configuration,  depending  on  the  time  interval 
[tn,tn+ 1],  for  which  the  total  energy  is  exactly  conserved.  The  implementation  of  this 
result  involves  the  actual  determination  of  the  generalized  mid-point  configuration  as  part 
of  the  algorithm.  We  show  that  this  can  be  easily  accomplished  by  solving  an  additional 
single  scalar  equation. 

iii.  Algorithms  incorporating  numerical  dissipation.  For  structural  dynamics,  high 
frequency  algorithmic  dissipation  is  often  regarded  as  a  desirable  property  of  a  time¬ 
stepping  algorithm.  A  number  of  modifications  of  the  classical  Newmark  family  of  al¬ 
gorithms  have  been  proposed,  typically  in  the  form  of  linear  multi-step  methods,  which 
introduce  high  order  dissipation  and  preserve  second  order  accuracy;  see  e.g.,  BATHE  & 
Wilson  [1973],  Park  [1975],  Wood,  Bossak  &  Zienkiewicz  [1981]  and  Bazzi  &  An- 
DERHEGGEN  [1982].  A  recent  review  is  contained  in  HOFF  &  Pahl  [1988a, b].  We  have 
presented  a  modification  of  the  o-method  of  HlLBER,  HUGHES  &  TAYLOR  [1977]  which 
involves  only  a  trivial  modification  of  the  preceding  Newmark-like  algorithms  for  nonlin¬ 
ear  shells.  In  fact,  the  proposed  extension  of  the  a-method  modifies  the  conservation 
form  of  the  momentum  equations  with  initial  conditions  at  the  beginning  of  the  time  step 
which  does  not  affect  the  linearization  of  the  algorithm.  We  showed,  however,  that  any 
departure  from  the  mid-point  rule  designed  to  introduce  numerical  dissipation  destroys 
the  conservation  properties  of  the  algorithm. 

It  is  well-known  (see  e.g.,  GOUDREAU  &  TAYLOR  [1973]),  that  the  mid-point  rule 
exhibits  no  numerical  dissipation  for  any  frequency  range.  In  fact,  for  the  linear  prob¬ 
lem,  the  algorithm  possess  two  complex  conjugate  roots  with  unit  modulus  (and  a  zero 
spurious  root)  which  bifurcate  into  a  double  real  root  of  unit  modulus  at  infinite  sam¬ 
pling  frequencies.  According  to  a  standard  result;  see  e.g.,  GEAR  [1971],  the  presence  of  a 
double  root  with  unit  modulus  leads  to  a  weak  instability  which,  for  the  mid-point  rule, 
can  occur  only  for  infinite  sampling  frequencies.  As  noted  by  CARDONA  &  GERARDIN 
[1979],  infinite  frequencies  can  occur  in  nonlinear  dynamics  if  kinematic  constraints  are 
enforced  via  Lagrange  multipliers,  since  they  have  associated  massless  degrees  of  freedom. 
It  appears,  however,  that  this  situation  is  not  relevant  to  the  problem  at  hand  since  the 
kinematic  constraint  which  requires  the  director  field  to  be  unit  is  enforced  exactly,  leading 
to  a  finite  element  method  with  five  degrees  of  freedom  per  node.  A  similar  observation 
applies  to  the  dynamics  of  nonlinear  rods  (see  SlMO  &  DOBLARE  [1990]).  High-frequency 
numerical  dissipation,  however,  does  result  in  added  stability,  since  the  spectral  radius  of 
the  algorithm  becomes  strictly  less  that  one  at  infinite  sampling  frequencies,  thus  prevent¬ 
ing  the  appearance  of  weak  instabilities  associated  with  the  mid-point  rule.  This  feature 
is  confirmed  by  our  numerical  simulations. 
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iv.  Finite  element  spatial  discretization.  Arbitrary  finite  element  interpolations  of 
the  semi-discrete  shell  equations  do’  not  necessarily  preserve  the  conservation  properties 
of  the  time-stepping  algorithm.  We  have  proved  that  for  nonlinear  shells  the  conservation 
properties  of  the  algorithm  are  preserved  by  the  spatial  discretization  if  the  following 
interpolations  are  adopted: 

iv.l.  Both  the  unit  director  constraint  and  the  linearized  director  constraint  condi¬ 
tions  are  enforced  exactly  at  the  nodal  points  of  the  finite  element  discretization. 
By  design,  these  conditions  are  satisfied  by  the  proposed  class  of  time-stepping 
algorithms  and  the  configuration  update  procedure  which  makes  use  of  the 
exponential  map. 

iv.2.  The  discrete  configuration  space  and  the  finite  dimensional  subspace  of  admissi¬ 
ble  variations  are  constructed  using  isoparametric  finite  element  interpolations. 

iv.3.  A  row  sum  approximation  is  used  only  for  the  inertia  terms  associated  with  the 
angular  velocity  of  the  director  field. 

Condition  iv.3  ensures  that  the  S1  invariance  property  of  the  shell  equations  with  respect 
to  drill  rotations  is  also  present  in  the  discrete  model  and  makes  it  possible  to  reduce  to  a 
five  degree  of  freedom  per  node  problem. 


3.  Concluding  Remarks 

We  believe  that  this  research  project  has  led  to  important  accomplishments  in  the 
mathematical  modeling  and  numerical  simulation  of  flexible  structures.  In  particular,  the 
aforementioned  developments  have  led  to  very  attractive  computational  methodologies 
for  the  large-scale  and  long-term  simulation  of  such  structures  undergoing  large  overall 
motions.  We  would  like  to  point  out  that  important  improvements  of  these  results  have 
been  obtained  in  the  development  of  conserving  time-stepping  algorithms  in  the  subsequent 
research  efforts  under  the  continued  support  of  the  AFOSE.  See  SlMO  &  Tarnow  [1992], 
Simo  &  Tarnow  [1994],  Simo,  Tarnow  &  Doblare  [1994],  Lewis  &  Simo  [1994]  in 
this  respect. 
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